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MATRIX METHOD OF DETERMINING THE LONGITUDINAL-STABILITY COEFFICIENTS AND 
FREQUENCY RESPONSE OF AN AIRCRAFT FROM TRANSIENT FLIGHT DATA! 


By JaĮmes J. Donecan and Henry A, PEARSON 


SUMMARY 


„4 matrix method is presented for determining the longitudinal- 
stability coefficients and frequency response of an aircraft 
from arbitrary maneuvers. The method is devised so that it 
can be applied to time-histury measurements of combinations 
of such simple quantities as angle of attack, pitching velocity, 
load factor, elevator angle, and hinge moment to obtain the 
over-all coefficients. Although the method has been devised 
primarily for the evaluation of stability coefficients which are 
of primary interest in most aircraft loads and stability studies, 
it can be used also, with a simple additional computation, to 
determine the frequency-response characteristics. The entire 
procedure can be applied or extended to other problems which 
can be ezpressed by linear differential equations. 


INTRODUCTION 


The longitudinal characteristics of an aircraft are often 
related by a second-order linear differential equation in 
which the aireraft is assumed to have freedom in pitch and 
in vertical motion; changes in forward velocity are so small 
that they can be neglected. In the evaluation of tail loads, 
the coefficients of the differential equation and the elevator 
forcing function are generally assumed to be known and the 
response is to be determined. In the evaluation of gust 
problems the response and the coefficients are assumed to 
be known and the foreing function is to be determined. By 
analogy in stability and control work, it is desirable to 
determine the restoring-force and’ damping-force coefficients 
from known forcing functions and responses. In case the 
damping is small enough to obtain the rate of decay (or 
logarithmic decrement) and period from the oscillation, the 
required damping and restoring coefficients are easily com- 
puted, Models employed in rocket-powered and drop tests 
can be and usually are so ballasted that such well-defined 
oscillations are obtained; however, the longitudinal oscilla- 
tions of piloted airplanes ordinarily are nearly critically 
damped and this analysis procedure cannot be applied. In 
any case, additional data and analysis are required to 
evaluate the control-effectiveness coefficients. 

Appreciable work has been done recently in the field of 
determining the frequency-response characteristics of air- 
eraft in flight and evaluating the stability coefficients from 





the frequency-response data. In general, the methods for 
determining these relationships have been to impose actually 
prescribed motions such as unit steps, triangular pulses, or 
sinusoidal motions to the elevator by means of special equip- 
ment and then to measure the responses. The theoretical 
methods for reducing such data are usually tailored to fit 
the prescribed elevator motion. References 1 and 2 present 
methods of treating input and output data by Fourier 
analysis to determine the frequency response. Compared 
with the direct sine-wave input method of evaluating the 
frequency response, these methods require less special equip- 
ment and flight time at the expense of additional compnta- 
tion. For the practical application of the Fourier transform 
method, certain restrictions are placed on the nature of the 
input and the resultant output motions: the motions must 
start from a trimmed steady-state condition and, at the end 
of the transient period, must approach either the original 
or the new steady or quasi-steady trim conditions. 

In view of the complications and limitations of existing 
methods of flight evaluation of stability coefficients and fre- 
quency response, development of a simple and less restricted 
flight test and associated analysis was considered desir- 
able. A matrix method for evaluating the longitudinal- 
stability coefficients of an aircraft directly from the input 
and output time histories corresponding to arbitrary con- 
trol motions has been derived in the present report. The 
frequency response and some of the stability derivatives 
may be evaluated once these coefficients are known. Àl- 
though this method was derived to determine the second- 
order longitudinal response of an aircraft, it can be applied 
to other systems which can be approximated by second- 
order differential equations; extension of the method to 
higher-order linear systems is also possible. 


SYMBOLS 
And combinations of aerodynamic parameters (sec 
f table I) 
b wing span, feet 
b: tail span, feet 
chord, feet 
Cy hinge-moment coefficient _f 
7 Per, 


t Supersedes NACA TN 2370, “Matrix Method of Determining the Longitudinal-Stability Coefficients and Frequency Response of an Aircraft From Transient Flight Data” by James J. 


Donegan and Henry A. Pearson, 1951. 
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Ch rate of change of hinge-moment coefficient with 
elevator angle (0€}/08) 

Ct lift coefficient (Lla S) 

Cr pitching-moment coefficient of airplane without 
horizontal tail (A4b/¢S?) 

Cm, pitching-moment coefficient of isolated horizontal 
tail surface 

g acceleration due to gravity, feet per second per 
second 

H hinge moment 

hy airplane radius of gyration about pitching axis, feet 

K empirical constant denoting ratio of damping mo- 
ment of complete airplane to damping moment 
produced by tail 

L lift, pounds 

m airplane mass, slugs (Wig) 

M pitching moment of airplane, foot-pounds 

n airplane load factor 

q dynamic pressure, pounds per square foot 

l vs 

Ge 

S wing area, square feet 3 

S: horizontal-tail area, square feet 

i time, seconds 

yV true velocity, feet per second 

W airplane weight, pounds 

Ti length from center of gravity of airplane to aero- 
dynamic center of tail (negative for conventional 
airplanes), feet 

Kı, Ko, 

KiB, 

Ks, Ka, dimensional constants occurring in equations 

K, Ks, (see table I) 

ra, KP, 

KV, K8, 

K% 

a wing angle of attack, radians 

ey tail angle of attack, radians 

Y flight-path angle, radians 

6 angle of pitch (a+y) 

ô elevator deflection, radians 

; de 

€ downwash angle, radians Ta a) 

nt tail efficiency factor (ado) 

$ phase angle between ineremental load factor and 
elevator deflection, degrees 

p mass density of air, slugs per cubic foot 

T dummy variable of integration 

o elevator angular velocity, radians per second 


The notations & and 6, a and @, and so forth, denote single 
and double differentiations with respect to time. 


a bar over letter represents maximum value 
lal bars on sides of symbol represent absolute value 











Relative wind.” 


aC 





AL, 


Tangent to 
fight path. 





Flight path 





FIGURE 1.—Sign conventions employed. Positive directions shown. 


Matrix notation: | 
Ht | rectangular matrix 


1] square matrix 

{} column matrix 

KAN integrating matrix (see table II) 
||.A]| matrix defined by equation (24) 
all transpose of |{A|| 

Subscripts: 

} denotes row elements in matrix 

j denotes column elements in matrix 
t tail 


LONGITUDINAL EQUATIONS OF MOTION 
ELEVATOR MOTION 


In this section the usual longitudinal equations of motion 
following an elevator motion are derived in such a manner as 
to obtain expressions between some of the simple combina- 
tions of variables which are measurable in flight: namely, 
angle of attack and elevator angle, pitching angular velocity 
and elevator angle, or load factor and elevator angle. The 
usual assumptions of linearity, small angles, no loss in air 
speed during the maneuver, and no flexibility are implied. 

As in reference 3, the differential equations of motion of 
an airplane due to a given elevator deflection may be written 
as (see fig. 1 for es 


mi Ve AagS— (5 E) 118, A5=0 (1) 


a ae z . y sie ò Ly K 
S? 


d Cn 
SA Aé E = ng z- Ad—mK,?§ =0 (2) 
dé dé b, 
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By use of the definitions 


A@=Ay+Aa 
é=7+e (3) 
@=y14 


equations (1) and (2) are reducible to the following second- 
order differential equation giving the relation between angle 
of attack and elevator angle: 


é+K,at+K,Ae=K,A5+Ki 3 (4) 


where the A’s are the constants for a given set of conditions 
and are defined in table I. The coefficient K, represents an 
effective aerodynamic-damping coefficient; K, represents an 
effective aerodynamic restoring-force coefficient; K; and K, 
represent effective elevator-control power coefficients. 

An alternative form of equation (4), expressing the relation 
between angle of pitch and elevator angle, may be obtained 
by inserting relations of equation (3) into equation (4) and 
noting from equation (1) that 


y= Aidat A, Ad 
y= dat Ag aed 
`g s$ 
sy= ai faa dt+ Ae Abdt 
La 0) 
where A; and A, are combinations of aerodynamic param- 


eters defined in table I. The equation obtained after 
these preceding substitutions are made is 


y t 
GLK, 6+K, A0=K,A 4K, f NE (6) 
0 


where (see table I) 
K;= GHKid: + Kady 
K= A K+ AK; 


From the following definition for load-factor increment 














=; dC ae - 
An=— =e, w WIS Aa npa K Anat (7) 
it follows that 
WS dy 
Ac=> Ie An =R 
da f 
¿= Wis. iny 
aC, "A (8) 
da q 
OWS „o Aas 
=I TA? 
da 


g, | ev? {den S oe Siri aac 
me 2m | da kj" da, "* ke 





TABLE I.—DEFINITION OF CONSTANTS OCCURRING 
IN EQUATIONS 





Con- 
stant 


K pV aC, Sir? Koy +22 § 
1 on Re 
2m| da, kë aa 


Definition 




















a K Sef 
a yp m IS 











i | ove[ Cr, Sa, Cmi Se dOr, dti Ka 2p r2S2 
K; | Me + 7: a 5 
dé kp’ dë bikè da, dô 2 mk,? 
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V7 






































K| Ki- K: pos zt de, l 
Caz Oar V Vae 
K; Cy de dCz pS x, 
A) gues) t a a i ee a E 
By Rias Oa, da da 2m Vr: 
1 
KS| ~—K. 
3 Ün 3 
| dC, 
da 
0 Sic A a 
E| ys“ 
The AY § term in equation (8) was found to be small and 
1 


is omitted in the subsequent derivation. 


Substituting the results from equation (8) into equation 


(4) yields another form expressing the relation between 
measured load-factor increment and elevator angle as 


fit Kit K,An=K;,As+Kg3 (9) 
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where (see table I) K; and Kg are now different forms of the 
effective control power coefficients. 


HINGE MOMENT 


The coefficients K, to K; occurring in equations (4), (6), 
and (9) are those associated with the measured elevator- 
motion case. The use of the relation 


(=o aote Aa (10) 


gives the solution for Aé as 


A= sr. (C— Sar Aas) (11) 


Dò 


The increment in tail angle of attack to be substituted in _ 


equation (11) is given by 





= de dCyupS A de 
sa| aa (1-Fe-“Ga Sin fa) °F dat) 
(12) 
so that 
Oran de dOr p S L 

s= a E da 2 = a) 

oui Geel 

i (Set) R me) 


In order to shorten the subsequent derivation for the hinge- 
moment case, the term K,é in equation (4) and its counter- 
parts in equations (6) and (9) are omitted. This effect. is 
usually small; however, each individual case should be exam- 
ined to see whether the term warrants dropping. 

A substitution of the value of Aé given by equation (13) 
into equations (4), (6), and (9) gives the following three 
differential equations for the same combination of variables 
with C, and its integral replacing Aé: 


àt K à+ KA= KC, (14) 
6+ KY 6+ Ke s— KO, 42s oan Cd (15) 
+K +K An=KeC, (16) 


INTEGRAL FORM OF EQUATIONS 


Although equations (4), (6), (9) and (14), (15), (16) could 
be used to evaluate the effective K coefficients from flight 
measurements of Aa, Å, and An together with measurements 
of elevator angle, stick force, or hinge moment, it is seen that 
several differentiations of the measured data would be re- 
quired. Inasmuch as a numerical differentiation process is 








inherently more inaccurate than the corresponding integra- 
tion process, the preceding equations are changed and 
rearranged so that either Aa, 6, or An, which are to be the 
measured values, appear as separate quantities on one side 
of the equation and the operations on these quantities 
appear on the other side. In integral form the rearranged 
equations are 


t $ r t r 
Kf Re ut+K, | ["sadrdt—K, f [“ssarde— 
Q JO RR t wD DO 


K,["sadt=—e (17) 
0 
t Tt 
Ki0-+Ks{edt—Ke | Asdi—Ky | f asdrdt=—0 (18) 
` a N At 


t t 7 ft fr 
K f andit f f Andrdt—K, | Í Aédrdi— 
C o Jo Jodo 


[3 
K, i Asem, (19) 


Ke | sadit K [faa drdt- Kë K K E E 
0 JO 0 JO, 
(20) 


t T 
Ki0+Ks | odi-K? f Cat—Ke $ K E E EE 
i (21) 


K, n andit Ka | [[andrat—Ke f n C drdi=—An 


(22) 


In principle to solve any one of these equations for the A 
coefficients, it is only necessary to tabulate the recorded 
values of the two basic variables (for example, in equation 
(19) the values of An and Aé) at a number of points t, ta, ta 
and so forth along a given time history and to perform the 
indicated integrations from ¢=0 up to the time of the recorded 
value t; A number of simultaneous equations containing 
the unknown K’s result which are then solved. The number 
of equations can vary from a minimum, in which the number 
of ordinates is equal to the number of unknown K's, to the 
case where there are more equations than unknowns. When 
the number of ordinates equals the number of unknown K’s, 
the usual methods of solving simultaneous equations may be 
used to obtain the A’s; however, when there are more equa- 
tions than unknowns, a least-squares method is required to 
reduce the equations. Since the best average value of the 
K’s is obtained when many points along the time history are 
used, a least-squares procedure is generally preferable. 

Although the integration indicated in equations (17) to 
(22) ean actually be performed graphically from the time 
histories, it is deemed better to express the equations in 
matrix form in order to enable a complete numerical solution 
to be made, 
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MATRIX FORM OF EQUATIONS 


Since the derivation in matrix form for any one of equations (17) to (22) is the same as for any other equation, only 
equation (19), involving measured load factor and elevator angles, is used. In matrix form the system of simultaneous 
equations obtained from reading the time history of the load factor n against elevator angle ô in an arbitrary pull-up may be 


written 


fs T [£] 
andt f andrat— f" f 
wd JO 


h d- f“ f asaz {z= 


i 








In shorter form this expression may be rewritten as 


|All 


¢ ` 
=į{—ån;} 


where the matrix |jA|| is in general a rectangular matrix; 
that is, for every time f, one equation or one row of the 
matrix |[A|| is obtained. The individual elements of matrix 
Al] are evaluated from the known values of incremental 
load factor and incremental elevator angle. As mentioned 
previously, the integration may be performed graphically 
but in the present case, use is made of the integrating 
matrices derived in reference 4. Thus, any element in the 


(>t 
rectangular matrix (equation (23)) such as | Andt or 
OD 


é r 
| f Andr di may be expressed in matrix form as follows: 
O t0 


jf, and} =C {ane} (25) 


f f andra =|| alf f an dil 


The integrating matrix ||C;|| as derived in reference 4 is 
given in table IT, with a time interval Af=0.1 second. It 
should be noted that a sufficient number of time intervals 
within the natural period being computed must be chosen to 
give a solution; usually the shorter the time interval chosen 
for the integrating matrix, the more accurate will be the 
final solution. 


i hie ande f” fanar di— N ‘ Abdrdt— T aadel 


fana f” f andrai— f? f asar di— f” adel | 


"Abd at— f Abdt 





—An, 
K, 
— Ån, 
K; 
— An. 
; “Adi K, na 








After the elements of the matrix A (equations (23) and 
(24)) have been determined either by applying the inte- 
grating matrix or by graphical integration, the method of 
least squares is applied to the solution of the system of 
simultaneous equations. In matrix notation the least- 
squares solution involves multiplication of matrix < by its 
transpose A’ so that equation (24) becomes 


[d'A] (Ki}={—A’ Any} (26) 
where the matrix [4’A] would be a 4 by 4 matrix for equa- 
tions (18) and (19). Equation (26) can now be arranged to 
be solved directly for the K’s by multiplying by the inverse 
matrix [AA] so that finally 


Alternately the system of simultaneous equations repre- 
sented by equation (26) can be solved for the values of K by 
any of the well-known methods of solving sets of simul- 
taneous equations, that is, by eliminating the variables or by 
using Crout’s method (reference 5). The derivation in 
matrix form of any of the other equations from (17) to (22) 
is similar to the plan given for equation (19) and, therefore, 
is not given. 


FREQUENCY RESPONSE 
As first derived by Cornell Aeronautical Laboratory 





(reference 6), the frequency response was measured by 
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TABLE IT.—INTEGRATING MATRIX [C1] 


{Based on 0.1-see intervals} 

































































t 0 G1 0.2 0.3 0.4 0.5 0.6 0.7 = 
0 0 0 0 0 0 0 0 0 0 7 
a - 042667 - 066667 | —. 008333 0 0 0 0 0 Q 
2 « 033333 + 133333 033333 0 0 0 0 0 0 = à 
3 + 033333 . 133333 076000 - 066667 | -—. 008333 0 nE MS.) 0 Q 
A . 083334 < 133333 066667 . 133333 033333 0 0 0 3 
5 « 033333 133333 066667 + 133333 075000 066667 | —. 008333 0 0 
6 - 033333 133333 066667 _ + 133333 - 066667 . 133333 033333 0 0 = 
ef « 033333 133333 066667 . 133333 . 066667 « 133333 075000 - 066667 | —. 008333 oe 
8 . 033333 133333 066667 . 133333 06686 . 133333, 066667 . 133333 ot 
9 + 033333 133333 066667 . 133333 06666 . 133333 006667 + 133333 
1.0 + 033333 133333 066667 . 133333 060667 . 133333 « 066667 . 133333 ae 
LA + 033333 133333 066667 . 133383 066667 . 133333 . 066667 . 133333 
1.2 + 033333 133333 066667 . 133333 066667 « 133333 066667 . 133333 - 066067 eae 
13 . 033333 « 133333 066667 . 133333 066667 « 133333 066667 - 133333 
1.4 - 033333 « 133333 066667 . 133333 . 06666 « 133333, 066667 . 133333 
3 5 - 033333 . 133333 066667 . 133333 . 005667 + 183333 | 066667 . 133333 H au 
17 eMecnasant Il) (Seacecte: (lt Weveeeece ff cckeweccs t|) “Zelesvee, BSE eee |) candglewe EAS 5 " ` 





actually subjecting the airplane to sinusoidal elevator mo- 
tions of various frequencies by means of specially constructed 
apparatus. From these results the coefficients Ki, Ko, and 
so forth, which are significant in control and loads work, 
could be determined provided the equation of motion was 
assumed, 

In the present instance since the coefficients K, and Ky 
are determined directly from the equation of motion, the 
corresponding relations are given so that the frequency 
response, which is significant in the design of stable autopilot 
systems, can also be determined. 

When a sinusoidal elevator motion has been assumed, 
then equation (9), omitting the minor effects of Kô, becomes 


n+Kyit+ Kk An=K;é sin wt (28) 


where An is the load-factor increment and w is the angular 
velocity of the elevator. Since equation (28) is a linear 
equation with constant coefficients, the steady-state solutions 
are of the form 

n=ñ sin (wt-+¢) 


r=Nw cos (wt -+ e) 


ü=— nw sin (wt+¢) 


(29) 


By a substitution of these relations into equation (28) the 
following equation is obtained: 


ña sin (wt-+ ¢)+ Ayre cos (ott p) + 


Kñ sin (wt+¢)=K;6 sin wt (30) 


whieh may be rewritten as 


a(K:— oð sin (wt-+ ¢)+- Kye cos (wt-+ ¢)=K;6 sin wt (31) 
or 


B sin (ot +¢+6)=K;6 sin of (32) 
where 
B= K,b=7(K,— 0 + (Kyo)? (338) 
and 
= —Kyw 2 
—e=d=—tan-!? =o 
e=o=tan Know (34) 


From equation (33) the amplitude ratio of load factor to 





elevator angle is seen to be 


ia 
5 


= x = K: 
VEF HEK 





(35) 








and the phase angle at various frequencies is given by 
equation (34). 

In the present case the values of Ki, Ko, and Ky would 
have been derived from the flight measurements and the 
values of w would be assigned. 

For the measured hinge-moments case the values of K,° 
and K? would be used instead of K, and Ag, and so forth. 
The complete frequency-response relations and transfer 
functions including all derivatives and integral of 6 for 
equations (4), (6), and (9) are given in the appendix. 


DETERMINATION OF AERODYNAMIC DERIVATIVES 


The various K coefficients determined from the meas- 
ured values may be termed effective coefficients and include, 
to some extent, effects of some nonlinearities, elasticity 
and effects of other variables which are omitted in the 
usual analysis. In addition, as may be scen from table 
I, the K coefficients are combinations of various quantities 
involving known geometric qualities, the conditions of the 
problem as well as aerodynamic derivatives. The stability 
coefficients given in table I are expressed in a form suitable 
to loads work. In usual stability calculations, these co- 
efficients are generally expressed in a simpler form where 
the number of aerodynamic variables are reduced and, as 
a result, the coefficients are more easily approximated. 
dC, dOn de 4Cr, 
da’ da’ da da,’ 





A total of 10 aerodynamic variables 


dCr, 10m, IC, DO, 

di’ dë’ da,’ 05 
the coefficients of table I. Although all the aerodynamic 
derivatives cannot be determined directly from the four basic — 
coefficients (namely, A), K2, Ag, and Ay), engineering approxi- 
mations of the more significant derivatives can be obtained 
if values are assigned to either some of the more accurately 
known derivatives or to those factors having least influence 
on the problem. 





12, and K appear in the definitions of 
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The factors having the least influence on the problem are 


K, az, and the derivative 





AC ny : . 
a which, respectively, allow 


for the contribution of wing-fuselage damping, tail efficiency, 
and moment due to tail camber to which average values 
can be assigned. A representative value of A is 1.25. 
Representative values of 4; range from about 1.2 to 0.8 
with the higher limit applying to propeller-driven airplanes 
operating at low speed and full power and the lower limit 
applying at high speed with the propeller braking. An 
average value for jets or at zero thrust for propeller-driven 
dn, 

dë 
can be obtained from existing wind-tunnel data or by using 
theoretical methods; —0.5 is an average value for tail 
surfaces. 

Since, as may be seen from table I, K; is directly propor- 

Y 


Cr 
qs an effective value of this derivative can 





airplanes is about 0.9. A representative value of 


tional to 





be determined directly from the definition of Ay. 

In order to determine consistent values of the remaining 
dC, dEn de Cu dC, 
da’ da’ da da,’ da,’ 


and a, further values must be assigned to several of the 





significant aerodynamic derivatives 


remaining derivatives. The derivatives chosen would natu- 
rally be those for which values could be obtained from other 
sources with the greatest degree of accuracy. 


EXAMPLES 


In order to illustrate the foregoing method as well as the 
consistency of results obtained with different sets of instru- 
mentation, typical examples are given using data obtained 
from three flights (referred to as flight 1, flight 2, and flight 3) 
of a high-speed medium jet bomber. For flight 1, the method 
of a computation is obtained in sufficient detail to enable a 
reader not too familiar with the mathematical details to 
reproduce similar results. Flight 1 is further divided into 
case | where data for An and Aé are used and case II where 
data for 6 and Aé are used. References 7 and 8 may be 
consulted for introductory discussions of least-squares and 
matrix methods. 

Figure 2 shows the measured time histories of velocity, 
altitude, incremental elevator displacement, incremental load 
factor, and incremental pitching velocity obtained during a 
push-down pull-up maneuver. By means of the values from 
figure 2, increments in load factor and elevator angle at 
0.1-second intervals have been tabulated in columns 2 and 3 
of table IIT. The elements of the 4 matrix (equations (23) 
and (24)) ere given in columns 4 to 7 of table IIL. Each 
element in these columns has been determined by performing 
the indicated integrations on the results given in columns 
2 and 3. In this instance the integrations have been 
performed by use of the previously mentioned integrating 
matrix derived in reference 4. This method is particularly 
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Figure 2.—Time histories of velocity, altitude, incremental elevator displacement, incre- 
mental load factor (computed and measured), and incremental pitching velocity (computed 
and measured) for flight 1 at Mach number 0.40. 


suitable when automatic computing machines are available. 

The elements of matrix A (equation (23)) which are given 
in columns 4 to 7 of table IIT indicate that with the Af 
spacing used there are 23 equations involving the four 
unknown values of K. In order to obtain the least-squares 
solution of these equations, the transpose || A||” of matrix |[]| 
is required. The transpose matrix is obtained by inter- 
changing the rows and columns of matrix ||]|. 

The product of the 4-row, 23-column transpose matrix by 
the 23-row, 4-column original matrix yields the 4-row, 
4-column matrix in the coefficients of K;. The resulting four 
simultaneous equations are then solved by any of the well- 
known methods of solving sets of simultaneous equations. 

By performing the preceding operations, the following 
values of K were obtained from the data listed in table IH: 

K, K: K; By 


3. 314221 T. 339706 — 119. 553905 5. 819025 


In order to show how well these computed values of Æ 
represent the original data, they have been reinserted into 
equation (19) along with the measured values of Aé to deter- 
mine calculated values of An. The computed curve is given 
by the dashed line in figure 2 of the plot of An against t. 
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i 2 3 : 5 6 | 7 
a 
Time. t Acceleration Elevator angle i ifr t fr i 
{sec} increment, increment, Aé Í An dt f f ân dr dt -Í { Aé dr dt -f Adi 
An (radians) a o 0 o o 

0 0 o g 0 
wt 054 . 046687 . 004050 . 000225 —. 000091 —. 002485 
2 —.054 . 077666 -005400 . 000720 —. 000622 —. 008814 
3 — 11 -082902 —. 002133 ! 000967 —. 001903 —. 016863 
4 —. 254 . 085084 —. 019667 | —. 000040 š —. 004008 —. 025292 
5 — 444 . 089010 —, 054950 —. 003632 —. 006969 —. 033990 
16 1588 . 093810 —. 106933 —. 011587 —. 010821 —. 043124 
T —. 784 -098610 —. 175658 —. 025559 —. 015610 —. 052741 
‘8 —.965 . 103846 —. 263233 —. 047347 —. 021386 —. 062860 
19 —1.122 . 106900 —, 367483 —. 078747 —. 028197 —. 073105 
1.0 1.291 - 109081 —. 488033 —. 121386 —. 036076 —. 084211 
1 —1. 462 . 108645 —. 626166 : =. 176978 —.045041 —.095004 
| 1.2 —1.575 . 108645 —. 778500 —. 247093 —. 055094 —. 105954 
i 1.3 1.704 -108645 —. 943233 —. 333111 —. 066233 —. 116830 
i 1.4 —1. 739 . 107336 —1. 116166 —, 436013 —. 078457 —. 127640 
1.5 —1.837 -053668 —1. 296083 —. 556599 —. 091666 —. 135681 
1.6 —1.801 . 003491 —1. 479099 —. 695333 —. 105416 —. 138490 
7 —1. 700 . 004800 —1, 654399 —. 852104 —. 119294 —. 138905 
1.8 —1. 569 —. 004799 —1. 818099 —1. 025826 | —. 133202 —. 139086 
1.9 —1.305 —.017017 —1. 961174 —1. 214978 —. 147083 —. 137970 
2.0 ~—1. 116 —, 026179 —2, 081599 —1. 417305 | —. 160760 —. 135785 
2.1 —. 369 —. 036851 —~2, 181332 —1. 630682 | —. 174186 —, 132600 
2.2 —. 364 —. 041887 —2, 253466 —1. 852652 j —. 187254 —. 128629 
2.3 —.315 —. 045814 —2, 296799 —2. 08034} —. 199905 =. 124284 

DISCUSSION AY argentine 


The same process as was used for the relations of An and Aé 
was also applied to the relations of 6 and Aé shown in figure 2. 
The tabular material corresponding to table III is not 
included; the values of K obtained, however, were as follows: 


Kı K, K; Ki 
3. 13167 8. 4123 —7. 6212 —12. 1967 


These values of K when reinserted into equation (18) 
resulted in the computed curve of Ô given by the dashed 
curve of figure 2. 

In addition to the preceding computations, several push- 
down pull-up maneuvers, made under similar conditions of 
altitude, weight, and center-of-gravity positions, were ana- 
lyzed to obtain the variation of several of the computed K’s 
with Mach number. In this analysis only, the measurements 
of An and Aé were used. The results obtained for three Mach 
numbers are shown in figure 3. The short parts of the curves 
shown are the expected variations in the K’s. Table I shows 
that A, should vary linearly with speed and the other values 
of K should vary parabolically. .The curves shown are 
merely guides adjusted to pass through zero and through the 
value of K at the 0.45 Mach number point. 

The values of K, and K; shown in figure 3 were also inserted 
into equations (34) and (35) to determine the corresponding 
curves of frequency response. The results are given in 
figure 4. 

In addition the values of K,, Ae, and frequency response for 
case I have been computed by using the definitions of table I 
and aerodynamic derivatives obtained from wind-tunnel 
tests. These results are also shown in figures 3 and 4. The 
aerodynamic derivatives were listed in an unpublished report 
by the North American Aviation, Inc. and were obtained in 
the Southern California Cooperative Wind Tunnel. 





If only the frequeney response is desired, it can be deter- 
mined without recourse to the equations of motion; however, 
if the stability coefficients are desired, it will be necessary to 
use the equations of motion as has been done in the present 
report. For either case several mathematical methods are 
available (references 1, 2, and 6) to obtain these required 
quantities and all methods, if carried far enough, should 
yield similar results. Thus, the present method is basically 
no more accurate than any other method; however, it has the 
advantage of simple instrumentation and experimental 
procedures but may require more extensive computation. 

As with other methods where linearity is a basic assump- 
tion, most consistent results are to be expected when the 
maneuvers are confined to the angle-of-attack region where 
linearity exists. In order for the outlined mathematical 
procedures to succeed, the maneuvers should cover as much 
of the linear range as possible in a short period of time and 
the portion of the maneuver considered should be confined 
to that portion where the integrals are increasing, This 
practice insures that the elements of the original matrix <4 
are all different and that the subsequent least-squares 
matrix [A’A] is not ill-behaved. Enough of the response 
time history should be taken to cover a good portion of the 
natural period of the system. A point worth noting in 
connection with the use of the equations is that zero time is 
assumed as being at the start of the maneuver when the 
airplane is in steady flight. Since the present method is not 
restricted by the final condition, it offers the possibility of 
performing an analysis on fragments of curves with the result 
that any variations in the constants. may be determined. 
In such an analysis two possibilities occur: (1) where the 
fragments considered start from a fixed initial condition and 
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FIGURE 3.— Variations of Ki and K: with Mach number. 


become successively longer, and (2) where the fragments are 
taken as consecutive. In the first case, the present method 
may be applied without any modification; in the seeond case, 
the equations must be altered to introduce the initial condi- 
tions for each fragment. These possibilities have not, 
however, been explored. 

In the derivation given herein, lag in downwash has been 
included (see equation (2)) but unsteady lift effects have not. 
References 9 and 10 show that for the present purposes the 
inaccuracy of omitting unsteady flow effects, exeept down- 
wash lag, is probably no greater than the inaccuracies in the 
original assumptions or of the experimental data. 

Other terms and other combinations of measurements 
might have been included in the derivations given—for 
instance, the equations are readily adapted to measurement 
of tail load and either airplane load factor, airplane angle of 
attack, or pitching angular velocity. Additional terms may 
have been ineluded to account for flexibility. Also it is 
possible, as for example in the case of the hinge-moment 
relations, to include additional terms to account for elevator 
moment-of-inertia effect, rate of elevator motion, and so 
forth in order to make the methods more inclusive. The 
inclusion of these further terms, however; generally requires 
additional X's to be evaluated and would only be justified 
when the assumptions implied in the basie equations of 
motion can be more closely approached and when the 
accuracy of measurements is high. Although the method 
had been applied herein to second-order differential equations, 
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FIGURE 4.—Airplane frequency response. 


it may be extended to higher-order equations with the limita- 
tion that too many integrations destroy the conditioning of 
the equations used in determining the coefficients (equation 
(26)) and make the equations difficult to work with. 

The results of the sample computations in which two 
different sets of instrumentation were used indicate an 
average difference between the respective K coefficients of 
about 10 pereent. The use of a least-squares method 
permits calculation of a probable error, which is an indication 
of how well the second-order system and the coefficients 
(computed on the basis of 0.1-second time intervals) fit the 
data. The expression used in computing the probable error 
is ` 
3E m 


vBu 





K:=0.6 í 15y ZE 
where Bi is the main diagonal term of [A'A], E is the 
difference between the computed and measured value of the 
variable, N is the number of cases considered in the least- 
squares procedure, and « is the number of variables deter- 
mined. This probable error has been calculated for case 
I and case IT and indicates an error of +0.3 in Ay and +0.5 
in Ay for the computations in which the accelerometer 
measurements were used. These values are contrasted with 
probable errors of £0.1 and +0.3 for the pitching-angular- 
velocity measurements. These probable errors are asso- 
ciated with the very small differences between the solid-line 





242 REPORT 1070- 
and dashed-line curves shown in figure 2. Greater accuracy 
may be obtained by increasing instrument accuracy, record- 
reading accuracy, and correcting original data for instrument 
errors. Further accuracy in the method may always be 
attained by using smaller time intervals. 

The results shown in figure 3 for the three flights investi- 
gated give some idea of the scatter to be expected between 
runs as well as the variation of the coefficients A, aud K; with 
Mach number. As might be expected from. the definition, 
Ky is seen to vary linearly with Mach number with little scat- 
ter. On the other hand, the values of K; either indicate a 
linear variation with Mach number or a scatter about the 
expected parabolic variation. 


The computed values of K, and K; (fig. 3) Spiros iroi 


the wind-tunnel data are in fair agreement with the flight- 
test values. For many engineering purposes this agreement 
may be adequate and probably typical of what might be ex- 
pected if wind-tunnel data were used at the design stage. 
Since all of the A’s are defined in table I, the dynamic longi- 
tudinal characteristies of an aircraft may be estimated in the 


-NATIONAL ADVISORY COMMITTEE FOR AERONAUTICS 


design stage by computing the A’s and inserting the values 
in the frequency-response relations given in the appendix, 


CONCLUDING REMARKS 


A matrix method has been presented for determining the 
longitudinal-stability coefficients and frequency response of 
an aircraft from an analysis of arbitrary maneuvers in which 
simple instrumentation is used. Errors in instrument aceu- 
racy and probable errors due to the use of a least-squares 
method are briefly discussed. Possible improvements in the 
method are discussed but, as of the present, it appears im- 
provements would be justified only for those cases where the 
basic assumptions are closely approached and where instru- 
ment accuracy is high. The method is equally applicable to 
other problems which can be expressed by second-order 
differential equations. 


LANGLEY AERONAUTICAL LABORATORY, 
NATIONAL Apvisory COMMITTEE FOR AERONAUTICS, 
LANGLEY Freup, Va., December 15, 1950. 





APPENDIX 
FREQUENCY-RESPONSE RELATIONS 


In the body of the report the phase angles and amplitude ratio were given for only the simplest case. 


The complete 


frequency-response relations and transfer functions for the equations involving ô and its derivatives are now presented. 
If D represents the differential operator d/dt, then the steady-state response due to a sinusoidal forcing function can be ob- 


tained by substituting iw for D in the transfer function. 





The following relations were developed by this procedure: 

















































































































Equation Transfer Function Phase Angle Amplitude Ratio 
wy Lo j = 7 A âg K; K,D ee 2 aE 1% -1 Kao Aa Ky? + Kyo? A 
at Kia i K Aa=kK,Ab+K,5 Ke OS K, | Gog =tan E atta K, re R, oft Kiu 
do ro) =tan! -Kw | Ac s Éo 
| a+ Kya+ K âa=K; Ab Aé DP?+K,D+K, ee Kı— o? Aé pe A E ae f 
ae ae - f',.,,40 KDK, {hee _, Kso x RETK 
j+ K0 KA0=K;44+ K; È Adia DDED +E) oo,=tan Ka +tan K: STR Koa a 
ae ea =K.. so pAn Kit ket Ky! -1 Kw -1 _ Kew a Ke HK Ko | 
nK n+ K An=K: át K4 Ks Be DI+K DLR; Pny = tan can e ERIT TORTENET 
“ i f An D > An Kot Ke a 
n+ Kynt+ K,An=K, A+ Ko | T lige on, =tan T! gttn a SIS Eun Kia 
nt+Kint+ K:ân=K,åô We DI KDIE, n= tan "Tee i VB oY + Kt 
| . 
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